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Abstract 

Background: Witli advances in next generation sequencing teclinologies and genomic capture tecliniques, 
exome sequencing lias become a cost-effective approacli for mutation detection in genetic diseases. However, 
computational prediction of copy number variants (CNVs) from exome sequence data is a challenging task. Whilst 
numerous programs are available, they have different sensitivities, and have low sensitivity to detect smaller CNVs 
(1-4 exons). Additionally, exonic CNV discovery using standard aCGH has limitations due to the low probe density 
over exonic regions. The goal of our study was to develop a protocol to detect exonic CNVs (including shorter CNVs 
that cover 1-4 exons), combining computational prediction algorithms and a high-resolution custom CGH array. 

Results: We used six published CNV prediction programs (ExomeCNV, CONTRA, ExomeCopy, ExomeDepth, 
CoNIFER, XHMM) and an in-house modification to ExomeCopy and ExomeDepth (ExCopyDepth) for computational 
CNV prediction on 30 exomes from the 1000 genomes project and 9 exomes from primary immunodeficiency 
patients. CNV predictions were tested using a custom CGH array designed to capture all exons (exaCGH). After this 
validation, we next evaluated the computational prediction of shorter CNVs. ExomeCopy and the in-house modified 
algorithm, ExCopyDepth, showed the highest capability in detecting shorter CNVs. Finally, the performance of each 
computational program was assessed by calculating the sensitivity and false positive rate. 

Conclusions: In this paper, we assessed the ability of 6 computational programs to predict CNVs, focussing on 
short (1-4 exon) CNVs. We also tested these predictions using a custom array targeting exons. Based on these 
results, we propose a protocol to identify and confirm shorter exonic CNVs combining computational prediction 
algorithms and custom aCGH experiments. 
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Background 

With advances in next generation sequencing technolo- 
gies and genomic capture techniques, exome sequencing 
has become a cost-effective approach for mutation de- 
tection in genetic diseases [1]. The availability of efficient 
and robust analysis tools such as Genome Analysis Toolkit 
(GATK) makes it possible to discover SNPs and indels 
using sequence data M^ith high sensitivity and specificity 
[2]. Even though CNV Copy Number Variant discovery in 
whole-genome sequence data is performed with high ac- 
curacy [3], a number of different algorithms with variable 
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specificities are available to detect CNVs in exome data [4]. 
Thus, selecting the correct algorithm or algorithm combin- 
ation has become a bottleneck of exome CNV prediction. 
Moreover, CNV detection algorithms have low specificity 
and sensitivity in predicting small CNVs (1 to 4 exons) [5]. 

While computational approaches have limitations in 
predicting CNVs in exome sequence data, standard array 
comparative genomic hybridization (aCGH) used for 
genome-wide high-resolution CNV detection also show 
restrictions in exonic CNV detection due to the low probe 
coverage over exonic regions. Thus, the detection of ex- 
onic CNVs using both computational and aCGH based 
methods remains a challenging task. 

The goal of our study was to develop a protocol to detect 
exonic CNVs (including CNVs that cover 1 to 4 exons) 
from exome sequencing data by combining computational 
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prediction algorithms and a high-resolution custom CGH 
array. In this study, we predicted CNVs in 30 exomes ob- 
tained from the 1000 genomes project [6] using six re- 
cently published CNV detection programs along with an 
in-house modified algorithm. Computational CNV predic- 
tions were then tested using a custom CGH array focused 
on exonic regions. Next, true CNVs were identified by 
comparing computational predictions with the results of 
the CGH array. With the experimental validation of the 
computational predictions, the sensitivity and false positive 
rate of each program, or program combination, was deter- 
mined. Results of the computational prediction demon- 
strated a wide range in both CNV count and size. Finally, 
we studied the clinical utility of the algorithms used in our 
study by computational prediction and array confirmation 
of CNVs in 9 exomes from primary immunodeficiency 
patients. 

Methods 

Data sets 

Alignment data (BAM files) required for computational 
CNV prediction were obtained from the public data re- 
pository of the 1000 genomes exome project which tar- 
gets the CCDS gene set in 2500 individuals [7]. All of 
the individuals are from the CEU population. DNA sam- 
ples of the same exomes used in computational CNV 
predictions were obtained from the Coriell Institute for 
Medical Research [8]. 

Exome sequencing 

Exome capture was performed on primary immunodefi- 
ciency patient DNA using the Agilent SureSelect Human 
All Exome capture kit v. 5 (Agilent Technologies), and 
captured libraries were sequenced on an lUumina HiSeq 
2000. Sequence alignment was performed with Novolign 
(V2.07.17) [9]. Next, using GATK (V2.4-9) [2], the initial 
BAM files were realigned and the base quality scores 
were recalibrated. After marking the duplicates with Pic- 
ard (VI. 74) [10], the final set of alignment data (BAM 
files) required for computational CNV prediction were 
generated. 

Computational CNV prediction 

The CNV analysis in our study has two main sections: 
computational CNV prediction and validation using a 
custom CGH array (Additional file 1: Figure SI). Com- 
putational CNV prediction was performed employing 6 
CNV prediction programs (Table 1) and a complete list- 
ing of the parameters used for all the programs is pre- 
sented in Additional file 1: Table SI. 

These programs use different statistical models in CNV 
calling. However, ExomeCopy and ExomeDepth use a 
similar statistical approach based on Hidden Markov Mo- 
dels. The main differences between these two programs 



Table 1 Exome CNV prediction programs used in the 
study 

Program Reference 

ExomeCNV https://secure.genome.ucla.edu/index.php/ExomeCNV_ 
User_Guide [11] 

CONTRA http://sourceforge.net/apps/mediawiki/contra-cnv/index. 
php?title=Main_Page [12] 

ExomeCopy http://www.bioconductor.org/pacl<ages/2.9/bioc/html/ 
exomeCopy.html [13] 

ExomeDepth http://cran.r-project.org/web/packages/ExomeDepth/ 
index.html [14] 

CoNIFER http://conifer.sourceforge.net/index.html [5] 

XHMM http://atgu.mgh.harvard.edu/xhmm/index.shtml [15] 

are the implementation of exon length normalization 
(SubdivideGRange method) in ExomeCopy [13] and gen- 
eration of aggregate reference set (select.reference.set 
method) in ExomeDepth [14]. As further discussed in re- 
sults and discussion, these two programs exhibited the 
lowest and highest stringency in CNV calling. Therefore, 
to reach improved stringency in CNV calling and to utilize 
methods implemented in both ExomeCopy and Exome- 
Depth, we designed an algorithm (ExCopyDepth) combi- 
ning these two programs (Flow and implementation of 
ExCopyDepth; Additional file 1: Figure S2). In ExCopy- 
Depth, ExomeCopy is used for exon length normalization 
(SubdivideGRange method in ExomeCopy) and to extract 
read count and GC% data for these normalized exonic re- 
gions (CountBamlnGRanges, GetGCcontent methods in 
ExomeCopy). Next, the most suitable reference dataset for 
each target exome was estimated using read count data 
from multiple samples (select.reference.set method in 
ExomeDepth). Finally, CNV regions were predicted using 
ExomeDepth (CallCNVs method). The initial tile length 
normalization performed in our method optimizes the 
generation of the aggregate reference set and minimizes 
the effect of GC content while CNV calling. 

Custom array design (exaCGH) 

To verify the ability to computationally call exonic QSTVs, 
a custom CGH array targeting exonic regions (exaCGH) 
was designed using the Agilent e Array web portal [16]. 
The design process was initiated by the Agilent High- 
Definition Probe (HD-Probe) search using a genomic in- 
terval list that specifies the 1000 genome exonic regions. 
As low probe scores are observed for GC-rich regions, it 
is difficult to design high-scoring probes for all exonic 
regions. In addition, Agilent guidelines recommend 150- 
200 bp minimum probe spacing for optimum probe 
hybridization. Therefore, in order to search high scoring 
probes and to reach the recommended spacing between 
probes, genomic intervals were designed by adding 800 bp 
flanking regions to either side of the exonic regions de- 
fined in the 1000 genomes exome project. Although the 
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incorporation of 800 bp flanking regions extends the 
probe search into intronic regions, this shows no negative 
efl'ect on detecting short CNVs since the probes are in- 
ternal to the CNV. A flowchart of the exaCGH array de- 
sign (data preparation, Agilent eArray HD-Search and 
evaluation of the search results) is presented in Additional 
file 1: Figure S3. The complete design, and capture regions 
of exaCGH, are available on request from the authors. 
Since probes covering the 1000 genomes exome do not 
occupy the whole IxlM array, the remaining features 
were filled by probes covering exonic, CDS (coding re- 
gions) and CCDS (conserved coding) regions specified in 
the GENCODE vl5 annotation [17]. Similar to the exome 
definition file used previously, 800 bp flanking regions 
were added to the 5' and 3' ends of the GENCODE v.l5 
regions (Additional file 1: Figure S3). After searching 
probes for the exon array, the microarray was generated 
according to eArray guidelines. 

We then compared exaCGH with two commercially 
available platforms used to detect CNVs: Agilent IxlM 
array and Aff)^metrix CytoScan HD. In order to compare 
these designs, we utilized exonic and CDS regions speci- 
fied in GENCODE V.15 and 800 bp flanking regions 
were added to the 5 ' and 3 ' end. Next, we evaluated the 
probe count distribution of each design (Additional file 1: 
Figure S4). Probe distributions clearly show that only ~ 
10% of probes of Affymetrix CytoScan HD and -1% of 
probes of Agilent IxlM mapped to GENCODE V.15 with 
at least 4 probes per exonic and CDS region (including 
1600 bp flanking regions). However, in the exaCGH de- 
sign, considering only regions containing at least 4 probes, 
83.9% of probes map to GENCODE v.15. The analysis of 
probe distribution demonstrated that exonic coverage 
(with at least 4 probes) of exaCGH is over 99.28 MB of 
GENCODEv.l5, which is significantly higher than the 
Agilent IxlM (3.5 MB) and Affymetrix CytoScan HD 
(6.67 MB). Moreover, to reduce the amount of noise in 
CNV detection, it is recommended to only call CNVs with 
at least 4 probes when analyzing aCGH results. Thus, the 
exaCGH array is a much improved design for detecting 
exonic CNVs due to the high probe density over exonic 
regions and the high exonic coverage. 

Next, exaCGH experiments were performed using 9 
DNA samples from the 1000 genomes exome project and 
9 DNA samples from primary immunodeficiency patients 
following Agilent protocol V. 6.3. Agilent Genomic Work- 
bench was used to call CNV regions which were detected 
by at least by four probes (with minimum average absolute 
log ratio for deletion and duplication > =0.25). 

Results and discussion 

Computational CNV prediction using published algorithms 

For computational CNV prediction we used six different 
programs and an algorithm (ExCopyDepth) combining 



ExomeCopy and ExomeDepth to detect CNVs from 30 
exomes. Since these prediction algorithms depend on a 
reference dataset, it is essential that exomes are captured 
and sequenced using the same technology. However, 
since the technologies used to generate the 30 exomes 
obtained from the 1000 genomes project can be grouped 
into three different categories (Additional file 1: Table 
S2), all the exomes were categorized accordingly and 
separate computational CNV predictions were performed 
for each group. 

Results of the computational algorithms show striking 
variation in the length and number of CNVs predicted 
by the different programs (Figure la, b). For example, 
CoNIFER has a count range of 1-248 and ExomeCopy 
has a count range of 158-1837. In order to directly com- 
pare programs at this level, we identified individual 
exons covered by CNVs (exCNVs) from each program. 

Analysis of exCNVs shows that the number of exCNVs 
predicted by computational programs differ considerably 
from each other (Figure 2a, b). For example, ExomeCopy 
predicted the highest mean number of exCNVs (1537) 
including 741 deleted and 795 duplicated exons. The 
lowest number of exCNVs was reported by ExomeDepth 
(average exonic CNV count 47, with 19 average number 
of deletions and 27 duplications). There are also large 
differences in the range of the number of exCNVs pre- 
dicted by each program (Figure 2a, b). When comparing 
reference datasets used in the programs (Additional file 1: 
Table SI), ExomeCNV is the only program that uses a sin- 
gle sample (an external exome, NA12282) as a reference. 
This makes ExomeCNV results heavily dependent on 
depth of coverage distribution of the reference exome 
used. Thus ExomeCNV results were not utilized for 
downstream analysis. 

ExCopyDepth 

As shown in previous studies, ExomeCopy demonstrated 
higher sensitivity compared to normalization and fast seg- 
mentation methods [18] in finding CNVs that overlap few 
exons [13]. A higher density over the 10 bp-10 kb range of 
CNV length distribution (Figure la) and the highest 
exCNV count compared to other programs (Figure 2a, b), 
may indicate the better performance of ExomeCopy in de- 
tecting shorter CNVs with the lowest stringency (which 
was tested and confirmed in the following stages of our 
study). Additionally, ExomeDepth has been optimized for 
detecting rare variants [14]. The lowest exCNV count 
observed for ExomeDepth compared to other programs 
(Figure 2a, b) may indicate the capability of ExomeDepth 
in detecting rare variants with higher stringency. 

Since we are interested in detecting shorter CNVs with 
improved stringency, we combined different methods 
implemented in each of these programs in a modified al- 
gorithm, aiming to reach higher stringency compared to 
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Figure 1 Count and length distributions of CNVs predicted by the programs in the study, (a) Length distribution of CNVs predicted by 
eacli program, (b) Count (number of CNVs) distribution of CNVs predicted by eacli program. 



ExomeCopy and lower stringency compared to Exome- 
Depth. We reasoned that the higher stringency (lower 
CNV count) compared to ExomeCopy would have lo- 
wer false positive CNV count compared to ExomeCopy, 
which was tested and confirmed in a later stage of our 
study. 

Evaluation of ExCopyDepth demonstrated an aver- 
age of 711 exCNVs per exome, while other programs 
showed an average exCNV counts between 47 and 
1537 (Figure 2a, b). The ExCopyDepth count (711) 
that falls away from the extreme ends of the exCNV 
count range (47-1537) exhibits the moderate strin- 
gency in CNV calling of ExCopyDepth. 

Identification of overlapping exCNVs 

exCNV counts (from 1000 genome exomes) were spread 
over a wide range (0-2500; Figure 2a, b). We reasoned 
that intersecting (overlapping) exCNVs predicted by 
multiple programs would have a narrower range, and 
thus be a better estimate of true CNVs, which appeared 
to be the case (Figure 2c, d). 

The maximum average number of overlapping de- 
letions and duplications (per program combination in 
Figure 2c, d) were 242 and 228 (from the intersection of 
ExomeCopy and ExCopyDepth), which were significantly 



lower than the highest average numbers resulting from 
individual programs (741 deletions and 795 duplications. 
Figure 2a, b). The number of overlapping exCNVs re- 
sulting from the intersection of three or more programs 
(combination 2 to 6 in Figure 2c, d) demonstrated a fur- 
ther decrease in exCNV count. The lower overlapping 
CNV counts from the intersection of ExomeCopy and 
ExCopyDepth compared to individual programs may 
imply lower false positive rates. We thus evaluated these 
predictions on a custom exon array. 

Validation of computational predictions using a custom 
exon array (exaCGH) 

Since standard aCGH platforms exhibit low exonic co- 
verage (Additional file 1: Figure S4), we designed a cus- 
tom array (exaCGH) with at least 4 probes per exon (see 
Methods). Agilent aCGH microarrays have previously 
been reported to have the highest sensitivity in detecting 
single-copy alterations between 1 and 49 kb [19]. We ran 
nine 1000 genomes samples (Additional file 1: Table S3) 
with exaCGH to vaUdate CNV predictions from each 
program. Additionally, since lower exCNV counts were 
shown by the intersection of ExomeCopy and ExCopy- 
Depth (Figure 2c, d), these results were also used for ex- 
perimental validation. 
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Figure 2 exCNVs predicted by programs and overlapping exCNVs predicted by program combinations. Each dot represents an individual 
exome. (a) Number of exonic duplications predicted by each program, (b) Number of exonic deletions predicted by each program, (c) Number 
of overlapping duplications predicted by each program combination, (d) Number of overlapping deletions predicted by each program 
combination. Program combinations (c and d): 1, ExomeCopy/ExCopyDepth; 2, ExomeCopy/ExCopyDepth/CONTRA; 3, ExomeCopy/ 
ExCopyDepth/CoNIFER; 4, ExomeCopy/ExCopyDepth/XHMM; 5, ExomeCopy/ExCopyDepth/ExomeDepth; 6, ExomeCopy/ExCopyDepth/ 
ExomeDepth/XHMM; 7, ExomeCopy/ExCopyDepth/CONTRA/XHMM; 8, ExomeCopy/ExCopyDepth/CoNIFER/XHMM. 
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For experimental validation, true positives (TPs) were 
defined as CNVs which were detected by computational 
methods and confirmed by exaCGH. Next, the average 
true positive CNV count (AvgTP; for nine 1000 genomes 
samples) per program was calculated (Table 2). Exome- 
Copy and ExCopyDepth have a higher AvgTP compared 
to the other algorithms (Table 2). When considering the 
low AvgTP we observed for CoNIFER, ExomeDepth and 
XHMM, it is important to note that these programs are 
optimized for the identification of rare variants in large 
exomes data sets. Therefore, the use of a relatively small 
sample collection may contribute to the reported low 
TP CNV count. 

False positive (FP) CNVs, defined as CNVs which were 
predicted only by the computational algorithms, were 
then identified. Next, these FP CNVs were compared to 
the exaCGH array design and further analyzed. In order 
to be conservative in FP identification, only CNV re- 
gions that have at least 4 probes in the exaCGH design 
were selected as the final set of FP CNVs. Then, the ratio 
between TP and FP CNVs (TP/FP ratio) was calculated 
for each program (Table 2). CoNIFER (0.92) presented 
the highest TP/FP ratio indicating the lowest FP count. 
ExomeDepth (0.67) and ExCopyDepth (0.34) showed sec- 
ond and third highest TP/FP ratio indicating lower FP 
CNV counts. The lowest TP/FP ratio (with highest FP 
counts) was exhibited by ExomeCopy and CONTRA 
(0.20). 

Next, we excluded CNV calls from the X and Y chro- 
mosomes and the TP/FP ratio was calculated for each 
program to study the effect of using male and female 
exomes in a single reference pool The resulting TP/FP 
ratios were higher than the previously calculated ratios 
and indicated lower FP CNV count for autosomal chro- 
mosomes. The difference in the TP/FP ratio demon- 
strated the effect of using male and female samples in 
the same reference pool Thus, splitting the reference 



exome collection based on gender (when large sample 
collections are available) will improve the CNV calling of 
X and Y chromosomes. 

When considering CONTRA results (Additional file 1: 
Table SI), it was observed that CONTRA predicts the 
coordinates of each exon with copy number alteration 
while all the other programs predict CNVs with one or 
multiple exons. Thus, in order to find CONTRA CNVs 
that contain one or multiple exons, adjacent exons were 
merged and considered as single CNVs, leaving isolated 
exons as CNVs with single exons. This set of CNVs are 
referred to as CONTRAmerged. Next, CONTRAmerged 
CNVs were compared to exaCGH results and TP, AvgTP 
and FP were identified. The TP CONTRAmerged CNV 
count was 67 (AvgTP: 7.5). CONTRAmerged CNV 
counts showed a clear decrease compared to the initial 
CONTRA results (AvgTP: 17.88). The FP CONTRA- 
merged count was 826 and the TP/FP ratio was 0.08. 
This is the lowest TP/FP ratio for all programs, and indi- 
cates that CONTRAmerged CNVs exhibits decreased per- 
formance in the implemented merge process. Since the 
next stages of our study are focused on further analysis of 
CNVs with one or multiple exons, CONTRA results were 
excluded from downstream analysis and CNVs predictions 
from other programs (without post-processing CNV calls) 
were utilized. 

Analysis of CNVs that cover 1-4 exons 

In order to evaluate the utility of computational pro- 
grams in detecting shorter CNVs (1-4 exons), we further 
analyzed the exon count of each predicted CNV. We 
counted the number of exons within each predicted 
CNV, and CNVs containing 1-4 exons were identified. 
The TP/FP ratio for each program was then calculated 
(Table 3). CoNIFER, ExomeDepth and XHMM did not 
identify any TP short CNVs. The intersection of Exome- 
Copy and ExCopyDepth showed the highest TP/FP ratio. 



Table 2 True positive (TP)/false positive (FP) CNV ratio predicted from each program 



Program 


AvgTP 


TP/FP ratio 


Total CNV count^ 


Average CNVs per sample^ 


CoNIFER 


1.33 (1.33)^ 


0.92 (1.09)2 


23 


2.56 


ExCopyDepth 


28.11 (15.86)^ 


0.34 (0.51)2 


422 


46.89 


ExomeCopy 


289.33 (226.0)^ 


0.20 (0.21)2 


11978 


1330.89 


ExomeDepth 


1.11 (0.78)^ 


0.67 (0.70)2 


17 


1.89 


Intersection of ExomeCopy / ExCopyDepth 


12.89 (7.33)^ 


0.28 (0.33)2 


218 


24.22 


XHMM 


4.44 (3.56)^ 


0.32 (0.53)2 


92 


10.22 


CONTRA 


17.88 (16.56)^ 


0.20 (0.20)2 


896 


99.56 



TP/FP ratio for each program was calculated using CNVs identified from 9, 1000 genomes samples run in both computational programs and exaCGH; Average true 
positive (AvgTP) = TP/9; TP/FP ratio = TP CNV count /FP CNV count; Average CNVs per sample = Total CNV count/9. 
^Average true positive calculated by excluding CNVs in X and Y Chromosomes. 
^TP/FP ratio calculated by excluding CNVs in X and Y Chromosomes. 

^Total number of CNVs predicted by each program excluding CNVs in X and Y Chromosomes (CNV counts for each program including X and Y Chromosomes are 
presented in Figure lb). 

"^Average CNVs per sample calculated from counts presented in total CNV count column in Table 2. 
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Table 3 True positive (TP), false positive (FP) and TP/FP ratio for short CNVs (1-4 exons) 



Program 


1000 genomes exomes 




Primary immunodeficiency patients 


TP 


FP 


TP/FP ratio 


TP 


FP 


TP/FP ratio 


CONIFER 


0 


0 




4 


6 


0.67 


XHMM 


0 


4 




7 


32 


0.22 


ExomeDepth 


0 


2 




50 


301 


0.17 


ExomeCopy 


914 


6591 


0.14 


161 


1076 


0.15 


ExCopyDepth 


63 


267 


0.24 


52 


669 


0.08 


Intersection of ExomeCopy/ExCopy Depth 


34 


76 


0.45 


44 


478 


0.09 



TP, True positive CNVs (CNVs identified both exaCGH and computational programs); FP, False positive CNVs (CNVs identified only by computational programs but 
not by exaCGH); TP/FP ratio. Ratio between true positive and false positive CNVs. 



demonstrating the ability of ExomeCopy and ExCopy- 
Depth to detect CNVs that cover 1-4 exons with the 
highest FP CNV count. 

Mutation identification using exomes from primary 
immunodeficiency patients 

Since we only tested nine 1000 genome samples with the 
exaCGH array, and we wanted to further test the utility of 
computational programs in detecting shorter CNVs, we 
then ran CNV prediction programs on nine exome sam- 
ples from patients with primary immunodeficiency dis- 
order (PID). Overlapping CNVs from ExomeCopy and 
ExCopyDepth were identified by performing an intersec- 
tion between ExomeCopy, ExCopyDepth results. These 
results were then mapped to the 1000 genomes exome 
definition file and CNVs which were covered by 1 to 4 
exons were identified. 

In parallel to computational CNV prediction, DNA 
samples from the patients were analyzed by exaCGH 
and CNV calling was performed using the same criteria 
applied to the 1000 genomes exaCGH experiments. The 
TP, FP and TP/FP ratio was calculated for CNVs that 
cover 1-4 exons for each program. As shown in Table 3, 
all the programs were able to detect short exonic CNVs. 
ExomeCopy (161), ExCopyDepth (52) and ExomeDepth 
(50) presented higher TP CNV counts compared to 
other two programs. Moreover CoNIFER showed the 
lowest FP count for detecting shorter CNVs for the PID 
exome collection. This study of shorter exonic CNVs 
using primary immunodeficiency patients exhibited the 
applicability of CoNIFER with lower FP CNVs, and Exo- 
meCopy, ExCopyDepth and ExomeDepth with higher 
TP CNV counts. 

The comparison of short CNV detection showed a 
higher performance for the intersection of ExomeCopy 
and ExCopyDepth in 1000 genomes samples compared 
to the PID group. We reasoned that sequence coverage 
differences between the two exome groups may contribute 
to the contrast in short CNV detection. The percentage of 
bases covered by at least 15 reads in 1000 genomes 



exomes ranges from 87%-93%, while PID exomes have a 
very narrow range (98%-99%). This higher variability of 
depth of coverage across the samples in 1000 genomes 
collection compared to PID group may contribute to the 
difference in the performance. In addition, CNV calling of 
1000 genomes exomes was performed by splitting the 30 
exomes into 3 groups (with 8, 9 and 13 exomes per group) 
when generating reference collections (Additional file 1: 
Table S2) whereas PID CNV calling was performed using 
a reference collection with 22 samples. When considering 
previous benchmark studies on ExomeDepth, CoNIFER 
and XHMM, these programs were tested using large refer- 
ence collections (eg. authors of CoNIFER have tested the 
program using 366 exomes [5]) and XHMM explicitly 
claims that it is designed to call CNVs from reference col- 
lections with at least 50 exomes [20]. Thus, in addition to 
the higher variation of read depth across the reference 
pool, we believe that the use of small sample collec- 
tions in 1000 genomes exomes also contribute to de- 
crease the performances of ExomeDepth, CoNIFER 
and XHMM compared to the intersection of Exome- 
Copy and ExCopyDepth. 

In order to further evaluate short CNV prediction of 
each program, we used CNVs identified from both ex- 
ome groups (1000 genomes exomes and PID patients) 
and examined the number of exons within each pre- 
dicted TP CNV. Here we compared the relative cumula- 
tive frequency distributions of exon count per TP CNVs 
predicted by computational algorithms. Figure 3 clearly 
shows that only a small proportion of the CNVs identi- 
fied by CoNIFER and XHMM contained 1-4 exons com- 
pared to ExomeCopy and ExCopyDepth. For example, - 
15% of CoNIFER and XHMM TP CNVs contain 1-4 
exons, while -40% of intersection of ExomeCopy and 
ExCopyDepth and -45% of ExCopyDepth and -60% of 
ExomeCopy identified TP CNVs with 1-4 exons. Analysis 
of the exon count distribution from each program clearly 
demonstrated the improved performances of ExomeCopy, 
ExCopyDepth and the intersection of ExomeCopy and 
ExCopyDepth in detecting shorter CNVs compared to 
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ExomeDepth, XHMM and CoNIFER. Hence, we further 
studied the statistical models and methods implemented 
in each program and detailed comparison of these algo- 
rithmic features affecting short CNV prediction is pre- 
sented in Additional file 1: Text SI. 



Analysis of the copy number state (CNS) predicted by 
each program 

We analyzed TP CNVs to study the consistency of com- 
putational algorithms in predicting the CNS of a genomic 
locus. To do this, we compared the computationally pre- 
dicted genomic regions from each exome (1000 genomes 
exomes and PID patient exomes) to the corresponding 



exaCGH result. Then we studied the consistency of pre- 
dicted CNS of each locus by intersecting the results of 
multiple program combinations (Table 4). The results 
show that ExomeCopy, ExCopyDepth, ExomeDepth and 
Conifer predicted the same CNS only for 50% of the 
regions detected by all 4 programs. For only 4 regions 
did all 5 programs agree on the CNS. In contrast, 91.2% 
of ExCopyDepth and ExomeDepth predictions and 
97.4% of ExomeCopy, ExCopyDepth and ExomeDepth 
predictions showed the same CNSs. The observed 8.8% 
and 2.6% difference in CNS is due to the difference in 
the segmentation implemented in ExCopyDepth and 
ExomeDepth. With the analysis of CNSs, it is important 
to note that computational algorithms were better at 



Table 4 Analysis of copy number state (CNS) predicted by different program combinations 



Program combinations 


Number of genomic regions 
with different CNS 


Number of genomic 
regions witii same CNS 


% of genomic regions 
withi same CNS 


ExomeCopy, ExCopyDepth 


0 


9 


100.00 


ExomeCopy, ExCopyDepth, ExomeDepth 


1 


37 


97.37 


ExomeCopy, ExCopyDepth, ExomeDepth, CoNIFER 


3 


3 


50.00 


ExomeCopy, ExCopyDepth, ExomeDepth, CoNIFER, XHMM 


0 


4 


100.00 


ExomeCopy, ExCopyDepth, ExomeDepth, XHMM 


1 


9 


90.00 


ExomeCopy, ExomeDepth 


0 


4 


100.00 


ExCopyDepth, ExomeDepth 


3 


31 


91.18 


ExCopyDepth, ExomeDepth, XHMM 


0 


4 


100.00 
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predicting copy number variable regions rather than 
actual CNSs. 

Evaluation of the performance of computational CNV 
prediction algorithms 

In order to assess the performance of the programs used, 
we calculated the sensitivity and false positive rate (FPR) 
using experimentally validated computational predictions. 
In addition to TP and FP CNVs, false negative (FN) CNVs 
were also identified from both exome groups for sensitiv- 
ity calculations. FN CNVs were considered as the CNVs 
which were detected only by the exaCGH array experi- 
ments but not predicted by computational programs. 

Next, sensitivity and FPR were calculated for each pro- 
gram (Figure 4). ExomeCopy had higher FPR compared 
to the other programs while CoNIFER had the lowest 
false positive rate. When considering the lower sensitiv- 
ities observed for ExomeDepth, CoNIFER and XHMM, 
it is important to note that these programs were opti- 
mized to detect rare CNVs. CoNIFER and XHMM were 
further optimized to detect CNVs containing 3 or more 
exons [15]. Optimizations implemented in these algo- 
rithms contribute to the high FN count and affects the 



sensitivity, as acknowledged in previous studies [4]. 
When comparing ExomeCopy and ExCopyDepth, sensi- 
tivity scores of ExCopyDepth were spread over a wider 
range (0.4-1). Overlapping CNVs identified from the 
intersection of ExomeCopy and ExCopyDepth indicated 
a slight decrease in the sensitivity and a clear decrease in 
the FP rate compared to ExomeCopy. For example, 7 
samples (38.9%) showed very high FPRs (FPR > 0.95) for 
ExomeCopy (Figure 4a) while only 1 sample (5.5%) 
showed FPRs over 0.95 for intersection of ExomeCopy 
and ExCopyDepth (Figure 4c). 

With the evaluation of the performance of CNV predic- 
tion algorithms, the clinical utility of these algorithms was 
tested in patients with primary immunodeficiency (PID). 
The PID patient exomes were studied for disease-related 
variants, and in one patient we identified homozygosity 
for a previously reported disease-causing single-nucleotide 
variant in the Fanconi-associated gene FANCA [21]. How- 
ever, there was no known family consanguinity and com- 
pound heterozygosity was expected. Further examination 
of the region surrounding this variant revealed a low 
relative read coverage and several succeeding homozy- 
gous variants, stretching from FANCA exon 26 to 37 




0.0 0.4 0.8 
False positive rate 




0.0 0.4 0.8 
False positive rate 




0.0 0.4 

False positive rate 




0.0 0.4 0.8 
False positive rate 




0.0 0.4 0.8 
False positive rate 




0.0 0.4 0.8 
False positive rate 



Figure 4 Sensitivity versus false positive rate for CNV prediction, (a) ExomeCopy. (b) ExCopyDepth. (c) Intersection of ExomeCopy and 
ExCopyDepth (overlapping CNVs predicted by ExomeCopy and ExCopydepth). (d) ExomeDepth. (e) CoNIFER. (f) XHMM. Sensitivity = true positive 
CNVs/(true positive CNVs + false negative CNVs). False positive rate = false positive CNVs/(false positive CNVs + true positive CNVs). 
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(NM_000135.2), suggesting a deletion at this locus. Ana- 
lysis of CNVs showed that ExomeCopy, ExCopyDepth, 
ExomeDepth and CoNIFER (all the programs except 
XHMM) predicted the deletion between exon 26-37 in 
FANCA, exaCGH validated the result, showing a deletion 
covering the same region in FANCA (Additional file 1: 
Figure S5), and confirmed the higher sensitivity observed 
for these programs. 

In summary, ExomeCopy, ExCopyDepth and the inter- 
section of ExomeCopy showed a higher performance in 
short CNV prediction compared to ExomeDepth, XHMM 
and CoNIFER (Figure 3). Moreover, all the computational 
prediction programs have high FPRs (Figure 4). Thus, in 
order to predict and to confirm short exonic CNVs, we 
propose a protocol with two stages combining computer 
predictions and custom CGH array (exaCGH) experi- 
ments. Stage one: initial CNV prediction by selecting the 
program or program combination depending on the target 
exome collection and user requirement. Stage two: va- 
lidation of computational predictions using exaCGH 
experiments. When selecting a program (or program 
combination) for short CNV detection, the intersection of 
ExomeCopy and ExCopyDepth should be used as it 
demonstrated better performance in detecting shorter 
CNVs compared to ExomeDepth, XHMM and CoNIFER 
(Figure 3) with lower FPR compared to ExomeCopy. Im- 
portantly, the intersection of ExomeCopy and ExCopy- 
Depth showed improved performance compared to other 
programs when relatively small sample collections (8-13 
exomes) were used for the reference pool or when read 
depth of the exomes in a reference pool varied over a wide 
range (TP/FP ratios of 1000 genomes exomes. Table 3). 
ExomeDepth can be selected to predict rare exonic CNVs 
with higher sensitivity compared to CoNIFER and XHMM 
(Figure 4d) when relatively large exome collections were 
available. However, it is important to note that short CNV 
detection of ExomeDepth is superior only to CoNIFER 
and XHMM (Figure 3). If the user requirement is to reach 
the lowest FPR in predicting CNVs (regardless of the 
lowest short CNV detection), CoNIFER can be selected 
(Figure 4e). As discussed earlier, all the programs demon- 
strated high FPR and showed inconsistency in predicting 
copy number state (Table 4). Thus, computationally pre- 
dicted CNVs can be validated using exaCGH experiments 
(stage two of our protocol) and true positive CNVs can be 
confirmed. 

In an era where many exome sequencing projects are 
moving toward detecting disease associated variants in 
lists of candidate genes, our proposed two stage protocol 
can be implemented to detect and confirm short CNVs. 
Stage one can be implemented to predict short CNVs in 
candidate genes, this overcomes the large FP rate when 
looking genome-wide. Stage two can then be used to val- 
idate predicted CNVs by exaCGH. 



Conclusions 

In summary, computational methods for CNV identifica- 
tion show clear variation in the number and size of pre- 
dicted CNVs. Evaluating the sensitivity and false positive 
rate of computational programs identified the highest 
sensitivity for ExomeDepth in detecting rare CNVs and 
the lowest false positive rate for CoNIFER. Following the 
analysis of short CNVs (1-4 exons), we designed a proto- 
col to identify shorter exonic CNVs by combining compu- 
tational prediction and custom aCGH methods, which has 
high relevance for mutation detection in candidate gene 
studies. 

Additional file 



Additional file 1: Figure SI. CNV discovery pipeline. Figure S2. Flow 
and implementation of the in-house designed ExCopyDepth algorithm. 
Figure S3. Agilent HD-Probe design and quality control. Figure S4. 
Probe count distribution of exaCGH, Agilent IxlM array and Affymetrix 
CytoScan HD. Figure S5. Identification and experimentally validation of a 
CNV in a disease-causing gene using computational predictions and CGH 
array. Table SI. Parameters used for each program. Table S2. Summary of 
the technologies used to derive the exomes used in the study. Table S3. 
1000 genomes sample IDs of the exomes run in our custom CGH array 
(exaCGH). Text SI. Comparison of the statistical mode and the statistical 
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